DischargeRouting.f90 Source File
Perform discharge routing
Source Code
!! Perform discharge routing !|author: <a href="mailto:giovanni.ravazzani@polimi.it">Giovanni Ravazzani</a> ! license: <a href="http://www.gnu.org/licenses/">GPL</a> ! !### History ! ! current version 2.3 - 2nd December 2024 ! ! | version | date | comment | ! |----------|-------------|----------| ! | 1.0 | 07/Oct/2016 | Original code | ! | 1.1 | 08/Feb/2023 | DischargePointInit and DischargePointExport added | ! | 1.2 | 12/Feb/2023 | routing parameters are no more stored in stream network. No more need of HydroNetwork module | ! | 1.3 | 25/Jan/2024 | module name changed from SurfaceRouting to DischargeRouting | ! | 1.4 | 26/Jan/2024 | module adjusted to cope with new module Diversions | ! | 1.5 | 10/Apr/2024 | observed reservoir downstream discharge read from file | ! | 1.6 | 11/Apr/2024 | observed diverted discharge read from file | ! | 1.7 | 18/Apr/2024 | discharge routing parameters assigned for different subbasins | ! | 1.8 | 08/May/2024 | Qin written to output instead of Qout, to fix the bug of zero discharge in outlet cell | ! | 1.9 | 06/Jun/2024 | log message bug fixed when diversion is outside channel| ! | 2.0 | 14/Jun/2024 | back to Qout written to output file. Qout is defined in last downstream cell| ! | 2.1 | 03/Jul/2024 | Diverted discharge from diversions included in Qlateral of Muskingum | ! | 2.2 | 04/Sep/2024 | pools variable moved to Reservoirs module | ! | 2.3 | 02/Dec/2024 | excess of water retained in snow added to Qlateral of Muskingum | ! !### License ! license: GNU GPL <http://www.gnu.org/licenses/> ! !### Module Description ! Route discharge in river network. ! ! River network includes different elements: ! ! - hillslope reaches ! - channel reaches ! - reservoirs ! - diversions ! ! The user sets the elements to simulate and their properties ! with a configuration file like in the following example ! !``` !export-channel-grid = 0 ! !masks-number = 2 # number of masks to assign channel parameters (at least 1 for the base-mask must exist) ! ![reservoirs] ! file = ./conf/reservoirs.ini ! dt = 10 # if dt = 0 reservoirs are not solved ! dt-out = 3600 #optional, default = 0 ! ![diversions] ! file = ./conf/diversions.ini ! dt = 0 # 0 = suppresses diversion simulation. otherwise dt is set equals to discharge routing dt ! dt-out = 3600 ! ![discharge-in] ! scalar = 0.0 ! ![discharge-out] ! scalar = 0.0 ! ![discharge-lat] ! scalar = 0.0 ! ![base-mask] ! channel-initiation-method = area #area, ask, curvature ! channel-initiation-threshold = 4000000 # ask 600000 m^2 area 3500000 ! hillslope-width = 200 # cross section width (m) ! hillslope-alpha = 45. # slope of trapezoidal section side bank (degree) ! hillslope-ks = 2 #Strickler roughness coefficient m^1/3 s^-1 ! ![mask1] ! file = ./data/subbasin1.asc ! format = esri-ascii ! epsg = 32633 ! channel-initiation-method = area #area, ask, curvature ! channel-initiation-threshold = 400000 # ask 600000 m^2 area 3500000 ! hillslope-width = 200 # cross section width (m) ! hillslope-alpha = 45. # slope of trapezoidal section side bank (degree) ! hillslope-ks = 3 #Strickler roughness coefficient m^1/3 s^-1 ! !Table Start !Title: channel properties !Id: base-mask !#threshold indicates the basin area below which values of parameters are applied !Columns: [count] [threshold] [width] [alpha] [ks] !Units: [-] [m^2] [m] [deg] [m^1/3s^-1] ! 1 5000000 5 45 20 ! 2 10000000 7 45 25 ! 3 15000000 10 45 30 ! 4 20000000 20 45 35 ! 5 30000000 25 45 40 ! 6 100000000 35 45 45 ! 7 500000000 50 45 45 ! 8 1000000000 65 45 45 ! 9 2000000000 80 45 45 !Table End ! !Table Start !Title: channel properties !Id: mask1 !#threshold indicates the basin area below which values of parameters are applied !Columns: [count] [threshold] [width] [alpha] [ks] !Units: [-] [m^2] [m] [deg] [m^1/3s^-1] ! 1 5000000 5 45 15 ! 2 10000000 9 45 20 ! 3 15000000 20 45 25 ! 4 20000000 25 45 30 ! 5 30000000 30 45 35 ! 6 100000000 35 45 40 ! 7 500000000 40 45 40 ! 8 2000000000 45 45 45 !Table End ! !``` ! MODULE DischargeRouting ! Modules used: USE DataTypeSizes, ONLY : & ! Imported Parameters: float, & short USE Loglib, ONLY : & !Imported routines: Catch USE GridLib, ONLY : & !imported definitions: grid_real, & grid_integer, & !imported routines: NewGrid, & GridDestroy, & ExportGrid, & ESRI_ASCII USE GridOperations, ONLY : & !Imported routines: GridByIni, & CRSisEqual, & CellArea, & !Imported operators: ASSIGNMENT (=) USE IniLib, ONLY: & !Imported routines: IniOpen, & IniClose, & SectionIsPresent, & KeyIsPresent, & IniReadString, & IniReadInt,& IniReadReal, & !Imported type definitions: IniList USE StringManipulation, ONLY : & !Imported routines: ToString, & StringCompact, & StringToUpper USE Morphology, ONLY : & !Imported rituines: DownstreamCell, & CheckOutlet USE Reservoirs, ONLY : & !imported definition: Reservoir, & !imported routines: InitReservoirs, & GetReservoirById, & CountReservoirs, & OutReservoirsInit, & OutReservoirs, & !imported variables: pools USE Diversions, ONLY : & !imported definition: Diversion, & !Imported routines: InitDiversions, & OutDiversionsInit, & GetDiversionById, & OutDiversions, & !Imported variables: dtDiversion, & dtOutDiversion, & diversionChannels USE Chronos, ONLY : & !imported definition: DateTime, & !Imported operators: OPERATOR (==), & OPERATOR (+), & ASSIGNMENT( = ) USE ObservationalNetworks, ONLY : & !Imported routines: ReadData, & ReadMetadata, & WriteData, & WriteMetadata, & CopyNetwork, & DestroyNetwork, & AssignDataFromGrid, & !Imported defined variable: ObservationalNetwork USE RoutingModels, ONLY : & !Imported routines: LevelPool, & DivertFlow, & MuskingumCungeTodini, & MuskingumCunge USE TableLib, ONLY: & !imported definitions: Table, & TableCollection, & !Imported routines: TableNew, & TableGetValue USE Irrigation, ONLY: & !imported variables: Qirrigation USE Snow, ONLY: & !imported variables: excessWaterSnowFlux USE Utilities, ONLY : & !Imported routines: GetUnit USE MorphologicalProperties, ONLY : & !imported variables: flowAccumulation, & flowDirection, & dem, & streamNetwork, & streamNetworkCreated USE RiverDrainage, ONLY : & !imported routines: ChannelInitiation IMPLICIT NONE !Public routines PUBLIC :: DischargeRoute PUBLIC :: InitDischargeRouting PUBLIC :: DischargePointInit !public declarations: TYPE(grid_real) :: Qin !! inlet discharge at time t (current) TYPE(grid_real) :: Qout !! outlet discharge at time t (current) TYPE(grid_real) :: Pin !! inlet discharge at time t-1 (previous) TYPE(grid_real) :: Pout !! outlet discharge at time t-1 (previous) TYPE(grid_real) :: Plat !! lateral discharge at time t-1 (previous) TYPE(grid_real) :: waterDepth !! water depth from thalweg (m) TYPE(grid_real) :: topWidth !! channel top width (m) INTEGER (KIND = short) :: dtDischargeRouting !Local routines PRIVATE :: DischargePointExport !Local declarations: INTEGER (KIND = short), PRIVATE, & PARAMETER :: MISSING_DEF_INT = -9999 REAL (KIND = float) , PRIVATE, & PARAMETER :: MISSING_DEF_REAL = -9999.9 INTEGER (KIND = short), PRIVATE, & PARAMETER :: SLOPES = 0 INTEGER (KIND = short), PRIVATE, & PARAMETER :: MCT = 1 !Muskingum-Cunge-Todini INTEGER (KIND = short), PRIVATE, & PARAMETER :: MUSKINGUM = 2 !Muskingum INTEGER (KIND = short), PRIVATE, & PARAMETER :: INSTANTANEOUS = 3 !Infinite celerity (Lake) INTEGER (KIND = short), PRIVATE, & PARAMETER :: MC = 4 !Muskingum-Cunge TYPE (ObservationalNetwork), PRIVATE :: sites TYPE (ObservationalNetwork), PRIVATE :: sitesDischarge INTEGER (KIND = short), PRIVATE :: fileUnitPointDischarge TYPE (DateTime), PRIVATE :: timePointExport TYPE (DateTime), PRIVATE :: timePoolsExport TYPE (DateTime), PRIVATE :: timeDiversionsExport INTEGER (KIND = short) :: dtOutReservoirs TYPE (grid_integer), PRIVATE :: channel TYPE (grid_integer), PRIVATE :: dams !!location of reservoirs in raster map TYPE (grid_integer), PRIVATE :: divChannels !!location of diversions in raster map TYPE (grid_real), PRIVATE :: manning !! Manning roughness (s / m^1/3) TYPE (grid_real), PRIVATE :: bankSlope !!river bank slope (degree) TYPE (grid_real), PRIVATE :: sectionWidth !!river section width (m) TYPE (TableCollection), PRIVATE :: channelParameters !======= CONTAINS !======= ! Define procedures contained in this module. !============================================================================== !| Description: ! Initialize surface routing SUBROUTINE InitDischargeRouting & ! (fileini, time, path_out, flowdirection, domain, & storage, netRainfall, dtrk ) IMPLICIT NONE !Arguments with intent in: CHARACTER (LEN = *), INTENT(IN) :: fileini !!name of configuration file TYPE(DateTime), INTENT(IN) :: time CHARACTER (LEN = *), INTENT(IN) :: path_out !!path of output folder TYPE(grid_integer), INTENT(IN) :: flowdirection TYPE(grid_integer), INTENT(IN) :: domain !Arguments with intent out: TYPE(grid_real), INTENT(OUT) :: storage !water stored in channel cell TYPE(grid_real), INTENT(OUT) :: netRainfall ! net rainfall rate [m/s] INTEGER(KIND = short), INTENT(OUT) :: dtrk !local declarations: TYPE (IniList) :: ini INTEGER (KIND = short) :: i, j, k !CHARACTER (LEN=5) :: ic !!option for initial condition TYPE(Reservoir), POINTER :: p !!point to current reservoir TYPE(Diversion), POINTER :: d !!point to current diversion channel CHARACTER (LEN = 30) :: channelInitiationMethod !!method to assign channel initiation REAL (KIND = float) :: channelInitiationThreshold !!channel initiation threshold (m2) REAL (KIND = float) :: hillslopelWidth !!hillslope section width (m) REAL (KIND = float) :: hillslopeKs !!hillslope roughness (m^1/3 s-1) REAL (KIND = float) :: hillslopeBankSlope !!hillslope bank slope (degree) REAL (KIND = float) :: area !! basin area (m2) REAL (KIND = float) :: channelValue INTEGER (KIND = short) :: exportChannelGrid !!set channel grid export REAL (KIND = float) :: scalar INTEGER (KIND = short) :: nmasks !! number of masks to assign routing parameters TYPE (grid_integer) :: maskTmp CHARACTER (LEN = 100) :: stringTmp !-------------------------------end of declaration----------------------------- !check stream network IF ( .NOT. streamNetworkCreated ) THEN CALL Catch ('error', 'DischargeRouting', & 'stream network is not available, check morphological properties' ) END IF !load configuration file CALL IniOpen (fileini, ini) !number of masks to assign discharge routing parameters IF ( KeyIsPresent ('masks-number', ini) ) THEN nmasks = IniReadInt ('masks-number', ini) ELSE CALL Catch ('error', 'DischargeRouting', & 'masks-number missing in configuration file' ) END IF !channel initiation !------------------ !allocate channel grid CALL NewGrid ( channel, domain, 0) !set option to export channel grid exportChannelGrid = IniReadInt ('export-channel-grid', ini) !first use base mask channelInitiationMethod = IniReadString ('channel-initiation-method', & ini, section = 'base-mask') channelInitiationThreshold = IniReadReal ('channel-initiation-threshold', & ini, section = 'base-mask') CALL ChannelInitiation (method = channelInitiationMethod, & threshold = channelInitiationThreshold, & mask = domain, flowAcc = flowAccumulation,& flowDir = flowDirection, dem = dem, & channel = channel ) ! modify channel initiation grid if additional masks are required DO k = 1, nmasks - 1 stringTmp = 'mask' // ToString (k) channelInitiationMethod = IniReadString ('channel-initiation-method', & ini, section = stringTmp) channelInitiationThreshold = IniReadReal ('channel-initiation-threshold', & ini, section = stringTmp) CALL GridByIni (ini, maskTmp, section = stringTmp) CALL ChannelInitiation (method = channelInitiationMethod, & threshold = channelInitiationThreshold, & mask = maskTmp, flowAcc = flowAccumulation,& flowDir = flowDirection, dem = dem, & channel = channel ) CALL GridDestroy (maskTmp) END DO !export channel grid IF ( exportChannelGrid == 1 ) THEN CALL ExportGrid (channel, 'channel.asc', ESRI_ASCII) END IF !create discharge parameter maps !------------------------------- !initialize empty maps CALL NewGrid (layer = manning, grid = domain, initial = 0. ) CALL NewGrid (layer = bankSlope, grid = domain, initial = 0. ) CALL NewGrid (layer = sectionWidth, grid = domain, initial = 0. ) !read channel parameter table(s) CALL TableNew ( fileini, channelParameters ) !assign parameters to base mask hillslopelWidth = IniReadReal ('hillslope-width', ini, section = 'base-mask') hillslopeKs = IniReadReal ('hillslope-ks', ini, section = 'base-mask') hillslopeBankSlope = IniReadReal ('hillslope-alpha', ini, section = 'base-mask') DO i = 1, domain % idim DO j = 1, domain % jdim IF ( domain % mat (i,j) /= domain % nodata ) THEN area = flowAccumulation % mat (i,j) * CellArea (domain,i,j) IF ( channel % mat (i,j) == 0 ) THEN !hillslope cell manning % mat (i,j) = 1. / hillslopeKs bankSlope % mat (i,j) = hillslopeBankSlope sectionWidth % mat (i,j) = hillslopelWidth ELSE !channel cell !set manning roughness coefficient CALL TableGetValue ( valueIn = area, & tables = channelParameters, & id = 'base-mask', keyIn = 'threshold', & keyOut ='ks', match = 'linear', & valueOut = channelValue, & bound = 'extendlinear' ) manning % mat (i,j) = 1. / channelValue !set river bank slope CALL TableGetValue ( valueIn = area, & tables = channelParameters, & id = 'base-mask', keyIn = 'threshold', & keyOut ='alpha', match = 'linear', & valueOut = channelValue, & bound = 'extendlinear' ) bankSlope % mat (i,j) = channelValue !set section width CALL TableGetValue ( valueIn = area, & tables = channelParameters, & id = 'base-mask', keyIn = 'threshold', & keyOut ='width', match = 'linear', & valueOut = channelValue, & bound = 'extendlinear' ) sectionWidth % mat (i,j) = channelValue END IF END IF END DO END DO !assign parameters to sub masks DO k = 1, nmasks - 1 stringTmp = 'mask' // ToString (k) hillslopelWidth = IniReadReal ('hillslope-width', ini, & section = stringTmp) hillslopeKs = IniReadReal ('hillslope-ks', ini, & section = stringTmp) hillslopeBankSlope = IniReadReal ('hillslope-alpha', ini, & section = stringTmp) CALL GridByIni (ini, maskTmp, section = stringTmp) DO i = 1, domain % idim DO j = 1, domain % jdim IF ( maskTmp % mat (i,j) /= maskTmp % nodata ) THEN area = flowAccumulation % mat (i,j) * CellArea (domain,i,j) IF ( channel % mat (i,j) == 0 ) THEN !hillslope cell manning % mat (i,j) = 1. / hillslopeKs bankSlope % mat (i,j) = hillslopeBankSlope sectionWidth % mat (i,j) = hillslopelWidth ELSE !channel cell !set manning roughness coefficient CALL TableGetValue ( valueIn = area, & tables = channelParameters, & id = stringTmp, keyIn = 'threshold', & keyOut ='ks', match = 'linear', & valueOut = channelValue, & bound = 'extendlinear' ) manning % mat (i,j) = 1. / channelValue !set river bank slope CALL TableGetValue ( valueIn = area, & tables = channelParameters, & id = stringTmp, keyIn = 'threshold', & keyOut ='alpha', match = 'linear', & valueOut = channelValue, & bound = 'extendlinear' ) bankSlope % mat (i,j) = channelValue !set section width CALL TableGetValue ( valueIn = area, & tables = channelParameters, & id = stringTmp, keyIn = 'threshold', & keyOut ='width', match = 'linear', & valueOut = channelValue, & bound = 'extendlinear' ) sectionWidth % mat (i,j) = channelValue END IF END IF END DO END DO CALL GridDestroy (maskTmp) END DO !initial condition !----------------- !input discharge IF (SectionIsPresent('discharge-in', ini )) THEN IF (KeyIsPresent ('scalar', ini, 'discharge-in') ) THEN scalar = IniReadReal ('scalar', ini, 'discharge-in') CALL NewGrid (Pin, domain, scalar) ELSE CALL GridByIni (ini, Pin, section = 'discharge-in') END IF ELSE !grid is optional: set to default = 0 CALL NewGrid ( Pin, domain, 0. ) END IF !output discharge IF (SectionIsPresent('discharge-out', ini )) THEN IF (KeyIsPresent ('scalar', ini, 'discharge-out') ) THEN scalar = IniReadReal ('scalar', ini, 'discharge-out') CALL NewGrid (Pout, domain, scalar) ELSE CALL GridByIni (ini, Pout, section = 'discharge-out') END IF ELSE !grid is optional: set to default = 0 CALL NewGrid ( Pout, domain, 0. ) END IF !lateral discharge IF (SectionIsPresent('discharge-lat', ini )) THEN IF (KeyIsPresent ('scalar', ini, 'discharge-lat') ) THEN scalar = IniReadReal ('scalar', ini, 'discharge-lat') CALL NewGrid (Plat, domain, scalar) ELSE CALL GridByIni (ini, Plat, section = 'discharge-lat') END IF ELSE !grid is optional: set to default = 0 CALL NewGrid ( Plat, domain, 0. ) END IF !allocate variables CALL NewGrid (layer = Qin, grid = domain, initial = 0. ) CALL NewGrid (layer = Qout, grid = domain, initial = 0. ) CALL NewGrid (layer = storage, grid = domain, initial = 0. ) !net rainfall grid CALL NewGrid (layer = netRainfall, grid = domain, initial = 0.) !water depth grid CALL NewGrid (layer = waterDepth, grid = domain, initial = 0.) !topwidth CALL NewGrid (layer = topWidth, grid = domain, initial = 0.) !diversions IF (SectionIsPresent('diversions', ini)) THEN dtDiversion = IniReadInt ('dt', ini, section = 'diversions') IF ( dtDiversion /= 0 ) THEN !set dtDiversion = dtDischargeRouting dtDiversion = dtDischargeRouting END IF IF (dtDiversion > 0) THEN CALL InitDiversions ( IniReadString('file', ini, & section = 'diversions') ) !create grid with diversion location CALL NewGrid (divChannels, domain, domain % nodata ) d => diversionChannels DO i = d % r j = d % c divChannels % mat (i,j) = d % id IF ( channel % mat (i,j) == 0 ) THEN CALL Catch ('warning', 'DischargeRouting', & 'diversion is not on channel: ', argument = ToString(d%id) ) END IF IF ( .NOT. ASSOCIATED (d % next) ) THEN !found last element of list EXIT END IF !set next diversion d => d % next END DO !initialize files for output IF ( KeyIsPresent ('dt-out', ini, 'diversions' ) ) THEN dtOutDiversion = IniReadInt ('dt-out', ini, 'diversions') IF ( dtOutDiversion > 0) THEN timeDiversionsExport = time CALL OutDiversionsInit ( path_out ) END IF END IF END IF ELSE !section is mandatory: stop the program if missing CALL Catch ('error', 'DischargeRouting', & 'error loading diversions: ' , & argument = 'section not defined in ini file' ) END IF !reservoirs IF (SectionIsPresent('reservoirs', ini)) THEN !create grid with reservoirs location CALL NewGrid (dams, domain, domain % nodata ) dtrk = IniReadInt ('dt', ini, section = 'reservoirs') IF (dtrk > 0) THEN CALL InitReservoirs (IniReadString('file', ini, & section = 'reservoirs'), & time, domain, pools) p => pools DO i = p % r j = p % c dams % mat (i,j) = p % id IF ( channel % mat (i,j) == 0 ) THEN CALL Catch ('warning', 'DischargeRouting', & 'reservoir is not on channel: ', argument = ToString(p%id) ) END IF IF ( .NOT. ASSOCIATED (p % next) ) THEN !found last element of list EXIT END IF p => p % next END DO !initialize files for output IF ( KeyIsPresent ('dt-out', ini, 'reservoirs' ) ) THEN dtOutReservoirs = IniReadInt ('dt-out', ini, 'reservoirs') IF ( dtOutReservoirs > 0) THEN timePoolsExport = time CALL OutReservoirsInit (pools, path_out) END IF END IF END IF ELSE !section is optional dtrk = 0 END IF !finished configuration, deallocate memory CALL IniClose (ini) END SUBROUTINE InitDischargeRouting !============================================================================== !| Description: ! route discharge in surface network SUBROUTINE DischargeRoute & ! (dt, time, flowdir, runoff, dtrk, riverToGroundwater, & groundwaterToRiver, storage ) IMPLICIT NONE !Arguments with intent(in): INTEGER(KIND = short), INTENT(IN) :: dt !!time step [s] TYPE(DateTime), INTENT(IN) :: time TYPE(grid_integer), INTENT(IN) :: flowdir !!flow direction TYPE(grid_real), INTENT(IN) :: runoff !! net rainfall [m/s] INTEGER(KIND = short), INTENT(IN) :: dtrk !!time step for reservoirs TYPE(grid_real), INTENT(IN) :: riverToGroundwater !!discharge !!from river to groundwater (m3/s) TYPE(grid_real), INTENT(IN) :: groundwaterToRiver !!discharge !!from groundwater to river (m3/s) !Arguments with intent inout TYPE(grid_real), INTENT(INOUT) :: storage !!volume stored on channel cell !local declarations: INTEGER(KIND = short) :: k, iin, jin, is, js REAL (KIND = float) :: ddx REAL (KIND = float) :: Qlateral !!lateral discharge (m3/s) REAL (KIND = float) :: QlateralChannel !!lateral discharge in a diversion channel (m3/s) REAL (KIND = float) :: Qdownstream !! discharge downstream a diversion channel (m3/s) REAL (KIND = float) :: Qdiverted !! discharge diverted from a diversion channel (m3/s) REAL (KIND = float) :: runoff_ij !!runoff of single cell REAL (KIND = float) :: wDepth REAL (KIND = float) :: tWidth REAL (KIND = float) :: area TYPE (Reservoir), POINTER :: p !!pointer to current reservoir TYPE (Diversion), POINTER :: d, dnest !!pointer to current diversion !-------------------------------end of declaration----------------------------- !reset Qin Qin = 0. !loop troughout hydro network DO k = 1, streamNetwork % nreach iin = streamNetwork % branch (k) % i0 jin = streamNetwork % branch (k) % j0 !follow the branch DO WHILE ( .NOT.( jin == streamNetwork % branch (k) % j1 .AND. & iin == streamNetwork % branch (k) % i1 ) ) !set runoff_ij runoff_ij = runoff % mat(iin,jin) !find downstream cell CALL DownstreamCell (iin, jin, flowdir % mat(iin,jin), & is,js, ddx, flowdir) ! looking for reservoir IF ( dtrk > 0) THEN IF ( dams % mat (iin,jin) /= dams % nodata ) THEN !set current reservoir p => GetReservoirById (list = pools, id = dams % mat (iin,jin) ) IF (time == p % tReadNewStage) THEN !read new stage value from file CALL ReadData (network = p % network, & fileunit = p % unit, time = time) !update target level p % stageTarget = p % network % obs (1) % value p % tReadNewStage = p % network % time !update time of next reading of target stage p % tReadNewStage = p % tReadNewStage + & p % network % timeIncrement END IF IF (p % dischargeDownstream .AND. & time == p % tReadNewDischargeDownstream ) THEN !read new downstream discharge value from file CALL ReadData (network = p % networkDischargeDownstream, & fileunit = p % unitDischargeDownstream, time = time) !update time of next reading of downstream discharge p % tReadNewDischargeDownstream = & p % tReadNewDischargeDownstream + & p % networkDischargeDownstream % timeIncrement END IF IF (p % dischargeDiverted .AND. & time == p % tReadNewDischargeDiverted ) THEN !read new diverted discharge value from file CALL ReadData (network = p % networkDischargeDiverted, & fileunit = p % unitDischargeDiverted, time = time) !update time of next reading of diverted discharge p % tReadNewDischargeDiverted = & p % tReadNewDischargeDiverted + & p % networkDischargeDiverted % timeIncrement END IF !reservoir routing Qin % mat(iin,jin) = Qin % mat(iin,jin) + & runoff_ij * & CellArea(runoff,iin,jin) CALL LevelPool (time, dt, dtrk, Pin % mat(iin,jin), & Qin % mat(iin,jin), p) !check and simulate diversion from reservoir IF ( p % bypassIsPresent ) THEN !divert flow from river d => p % bypass CALL DivertFlow (time, d, p % Qout, p % Qout ) !overwrite diversion inlet discharge when observation is available IF ( p % dischargeDiverted ) THEN IF ( p % networkDischargeDiverted % obs (1) % value /= & p % networkDischargeDiverted % nodata ) THEN d % QinChannel = p % networkDischargeDiverted % obs (1) % value END IF END IF !route discharge into channel Qlateral = 0. CALL MuskingumCungeTodini ( dt, & d % channelLenght, & d % QinChannel, & d % PinChannel, & d % PoutChannel, & Qlateral, Qlateral, & d % channelWidth, & d % channelBankSlope, & d % channelSlope, & d % channelManning, & d % QoutChannel, & wDepth, tWidth ) !copy data of current step for next step d % PinChannel = d % QinChannel d % PoutChannel = d % QoutChannel END IF !set Qout SELECT CASE ( p % typ ) CASE ( 'on' ) IF ( p % dischargeDownstream ) THEN IF ( p % networkDischargeDownstream % obs (1) % value /= & p % networkDischargeDownstream % nodata ) THEN !overwrite reservoir downstream discharge p % Qout = p % networkDischargeDownstream % obs (1) % value END IF END IF Qout % mat(iin,jin) = p % Qout CASE ( 'off' ) !off-stream reservoir DA RIVEDERE!! !compute off-stream pool outflow CALL TableGetValue ( valueIn = p % stage, tab = p % geometry, keyIn = 'h', & keyOut ='Qout', match = 'linear', & valueOut = p % Qout_off, bound = 'extendlinear' ) p % Pout_off = p % Qout_off Qout % mat(iin,jin) = p % Qout !assign outflow to pool out cell !IF (p % rout /= MISSING_DEF_INT .AND. p % cout /= MISSING_DEF_INT ) THEN !Qin % mat(p % rout,p % cout) = Qin % mat(p % rout,p % cout) + & ! p % Qout_off !END IF END SELECT Pin % mat(iin,jin) = Qin % mat(iin,jin) Pout % mat(iin,jin) = Qout % mat(iin,jin) jin = js iin = is CYCLE END IF END IF !Muskingum-Cunge-Todini area = CellArea(runoff,iin,jin) !set Qlateral Qlateral = runoff_ij * area !remove irrigation IF ( ALLOCATED (Qirrigation % mat) ) THEN Qlateral = Qlateral - Qirrigation % mat (iin, jin) END IF !include river-groundwater exchange IF ( ALLOCATED ( riverToGroundwater % mat ) ) THEN IF ( riverToGroundwater % mat (iin,jin) /= & riverToGroundwater % nodata ) THEN Qlateral = Qlateral + groundwaterToRiver % mat (iin,jin) - & riverToGroundwater % mat (iin,jin) END IF END IF !add excess of water retained in snow IF ( ALLOCATED ( excessWaterSnowFlux % mat ) ) THEN Qlateral = Qlateral + excessWaterSnowFlux % mat (iin, jin) * area END IF !add outflow from by-pass channel or off-stream basin p => pools DO WHILE (ASSOCIATED(p)) !loop trough all reservoirs IF ( p % typ == 'off' ) THEN IF ( iin == p % rout .AND. jin == p % cout ) THEN Qlateral = Qlateral + p % Pout_off END IF END IF IF (p % bypassIsPresent) THEN IF ( iin == p % bypass % rout .AND. jin == p % bypass % cout ) THEN Qlateral = Qlateral + p % bypass % PoutChannel END IF END IF p => p % next END DO !remove diverted discharge from diversion channel IF (dtDiversion > 0) THEN IF ( divChannels % mat (iin,jin) /= divChannels % nodata ) THEN !divert flow from river d => GetDiversionById (list = diversionChannels, & id = divChannels % mat (iin,jin) ) CALL DivertFlow (time, d, Qin % mat(iin,jin), Qdownstream ) Qdiverted = Qin % mat (iin, jin) - Qdownstream !update Qlateral Qlateral = Qlateral - Qdiverted !route discharge into channel QlateralChannel = 0. CALL MuskingumCungeTodini ( dtDiversion, & d % channelLenght, & d % QinChannel, & d % PinChannel, & d % PoutChannel, & QlateralChannel, QlateralChannel, & d % channelWidth, & d % channelBankSlope, & d % channelSlope, & d % channelManning, & d % QoutChannel, & wDepth, tWidth ) !copy data of current step for next step d % PinChannel = d % QinChannel d % PoutChannel = d % QoutChannel END IF END IF !add outflow from diversion channel d => diversionChannels DO WHILE (ASSOCIATED(d)) !loop trough all diversion channels IF ( iin == d % rout .AND. jin == d % cout ) THEN ! As the flood routing is solved from upstream to downstream ! PoutChannel contains the current or the previous time step discharge according to the ! Horton order of the outlet section respect to the horton order of onlet section Qlateral = Qlateral + d % PoutChannel END IF d => d % next END DO CALL MuskingumCungeTodini ( dt, ddx, Qin % mat(iin,jin), & Pin % mat(iin,jin), & Pout % mat(iin,jin), & Qlateral, Plat % mat(iin,jin), & sectionWidth % mat (iin,jin), & bankSlope % mat (iin,jin), & streamNetwork % branch (k) % slope, & manning % mat (iin,jin), & Qout % mat(iin,jin), & waterDepth % mat(iin,jin), & topWidth % mat(iin,jin) ) storage % mat (iin,jin) = storage % mat (iin,jin) + & ( Qin % mat(iin,jin) + Qlateral - Qout % mat(iin,jin) ) * & dt / CellArea(Qin,iin,jin) IF ( storage % mat (iin,jin) < 0. ) THEN !storage % mat (iin,jin) = 0. !Qin % mat(iin,jin) = 0. !Qout % mat(iin,jin) = 0. END IF IF (isnan (Qout % mat(iin,jin)) ) THEN Qout % mat(iin,jin) = Pout % mat(iin,jin) END IF !computed Qout becomes Qin for the downstream section. !Take account of junctions, sum to existing discharge Qin % mat(is,js) = Qin % mat(is,js) + Qout % mat(iin,jin) !store previous values for next time step Pin % mat(iin,jin) = Qin % mat(iin,jin) Pout % mat(iin,jin) = Qout % mat(iin,jin) Plat % mat(iin,jin) = Qlateral !store cell indexes and select downstream cell for next loop jin = js iin = is !check next cell is outlet cell CALL DownstreamCell (iin, jin, flowdir % mat(iin,jin), & is,js, ddx, flowdir) IF ( CheckOutlet (iin,jin,is,js,flowdir) ) THEN !next cell will not be computed, set approximate value of Qout Qout % mat (iin, jin) = Qin % mat (iin, jin) + & runoff % mat(iin,jin) * CellArea(runoff,iin,jin) END IF END DO END DO !write results IF ( time == timePointExport ) THEN CALL DischargePointExport ( time ) timePointExport = timePointExport + sitesDischarge % timeIncrement END IF IF ( time == timePoolsExport ) THEN CALL OutReservoirs ( pools, time, Qin, Qout ) timePoolsExport = timePoolsExport + dtOutReservoirs END IF IF ( time == timeDiversionsExport ) THEN CALL OutDiversions ( diversionChannels, time, Qin, Qout ) timeDiversionsExport = timeDiversionsExport + dtOutDiversion END IF RETURN END SUBROUTINE DischargeRoute !============================================================================== !| Description: ! Initialize export of point site data SUBROUTINE DischargePointInit & ! ( pointfile, path_out, time ) IMPLICIT NONE !Arguments with intent (in): CHARACTER (LEN = *), INTENT(IN) :: pointfile !!file containing coordinate of points CHARACTER (LEN = *), INTENT(IN) :: path_out !!path of output folder TYPE (DateTime), INTENT(IN) :: time !!start time !local declarations INTEGER (KIND = short) :: err_io INTEGER (KIND = short) :: fileunit !-------------------------end of declarations---------------------------------- timePointExport = time !open point file fileunit = GetUnit () OPEN ( unit = fileunit, file = pointfile(1:LEN_TRIM(pointfile)), & status='OLD', iostat = err_io ) IF ( err_io /= 0 ) THEN !file does not exist CALL Catch ('error', 'DischargeRouting', 'out point file does not exist') END IF !Read metadata CALL ReadMetadata (fileunit, sites) !check dt IF (.NOT. ( MOD ( sites % timeIncrement, dtDischargeRouting ) == 0 ) ) THEN CALL Catch ('error', 'DischargeRouting', & 'dt out sites must be multiple of dtDischargeRouting ') END IF CLOSE ( fileunit ) !create virtual network and initialize file for output fileUnitPointDischarge = GetUnit () OPEN ( unit = fileUnitPointDischarge, & file = TRIM(path_out) // 'point_discharge.fts' ) CALL CopyNetwork ( sites, sitesDischarge ) sitesDischarge % description = 'discharge data exported from FEST' sitesDischarge % unit = 'm3/s' sitesDischarge % offsetZ = 0. CALL WriteMetadata ( network = sitesDischarge, & fileunit = fileUnitPointDischarge ) CALL WriteData (sitesDischarge, fileUnitPointDischarge, .TRUE.) ! destroy sites CALL DestroyNetwork ( sites ) RETURN END SUBROUTINE DischargePointInit !============================================================================== !| Description: ! Export of point site data SUBROUTINE DischargePointExport & ! ( time ) IMPLICIT NONE !Arguments with intent(in): TYPE (DateTime), INTENT (IN) :: time !local declarations: INTEGER (KIND = short) :: i !-------------------------end of declarations---------------------------------- !set current time sitesDischarge % time = time !populate data !CALL AssignDataFromGrid (Qin, sitesDischarge ) !debug CALL AssignDataFromGrid (Qout, sitesDischarge ) !write data CALL WriteData (sitesDischarge, fileUnitPointDischarge) RETURN END SUBROUTINE DischargePointExport END MODULE DischargeRouting